function [ subjectinfo ] = analyzeVisualResponsivenessBanded_Calculate( subjectinfo )
%ANALYZEVISUALRESPONSIVENESSBANDED_CALCULATE Summary of this function goes here
%   Detailed explanation goes here

% Clear output window.
clc;
if ~isa(subjectinfo, 'SubjectInfo')
    error('Invalid subjectInfo.');
end
taskname = 'analyzeVisualResponsivenessBanded_Calculate';
% CheckTask(subjectinfo, taskname);
subjectinfo = DoWork(subjectinfo);
% CleanUp(subjectinfo, taskname);
end

%% Task related support
function CheckTask(subjectinfo, taskname)
% Check the task status
myTaskIndex = TaskInfo.FindTask(subjectinfo, taskname);
if (subjectinfo.Tasks(myTaskIndex).IsDone)
    disp(subjectinfo);
    disp(['Task "' taskname '" has already been completed.']);
    pause(5);
    return
end
if (subjectinfo.Tasks(myTaskIndex).HasDependancy)
    myDepTaskIndex = TaskInfo.FindTask(subjectinfo, subjectinfo.Tasks(myTaskIndex).Dependancy);
    if (subjectinfo.Tasks(myDepTaskIndex).IsToDo)
        error(['Task "' taskname '" depends on ' subjectinfo.Tasks(myDepTaskIndex).Taskname ', which hasn''t been completed yet.']);
    end
end
disp(['Running task "' taskname '"...']);
end


function CleanUp(subjectinfo, taskname)
myTaskIndex = TaskInfo.FindTask(subjectinfo, taskname);
% Set the task status
subjectinfo.Tasks(myTaskIndex) = subjectinfo.Tasks(myTaskIndex).SetComplete;
disp(['Task "' taskname '" completed.']);

% Save the subject information
subjectinfo.Save;
disp(subjectinfo);
pause(2); % Wait 2 seconds before continueing
end


%% Actual work..
function [subjectinfo] = DoWork(subjectinfo)

% Steal the epoch-work from Visual Responsiveness
subjectinfo.Analyses.VisualResponsivenessBanded = subjectinfo.Analyses.VisualResponsiveness;

% Configure the bands
data.bands = cell(1,7);
data.bands{1}.foi =   4:0.1:  7;    data.bands{1}.description = '4-7Hz (theta)';
data.bands{2}.foi =   7:0.1: 13;    data.bands{2}.description = '7-13Hz (alpha)';
% data.bands{3}.foi =   7:0.1: 15;    data.bands{3}.description = '7-15Hz (broad alpha)';
data.bands{3}.foi =  13:0.2: 30;    data.bands{3}.description = '13-30Hz (beta)';
% data.bands{5}.foi =  15:0.5: 40;    data.bands{5}.description = '15-40Hz (broad beta)';
data.bands{4}.foi =  30:0.5: 50;    data.bands{4}.description = '30-50Hz (low gamma)';
data.bands{5}.foi =  50:0.5: 90;    data.bands{5}.description = '50-90Hz (middle gamma)';
data.bands{6}.foi =  90:1.0:150;    data.bands{6}.description = '90-150Hz (high gamma)';
data.bands{7}.foi = 150:1.0:250;    data.bands{7}.description = '150-250Hz (very high gamma)';
timeStep = 0.05;
% timeStep = 1; 
% timeStep = 5;

% Process each set
% TODO: Make this a parfor!
for x = 1:length(subjectinfo.Analyses.VisualResponsivenessBanded.Sets)
    wavedataFilename = [subjectinfo.Analyses.VisualResponsivenessBanded.Filepath subjectinfo.PrimaryPrefix ' Wavedata for VisResp Banded Set ' num2str(x) '.mat'];
    
    [dataRest electrodeLabels] = LoadSet(subjectinfo.Analyses.VisualResponsivenessBanded.Sets(x).DataFilenameRest);
    dataStimulus = LoadSet(subjectinfo.Analyses.VisualResponsivenessBanded.Sets(x).DataFilenameStimulus);
   
    for b = 1:length(data.bands)
        waveRest = DoFreqAnalysis(dataRest, data.bands{b}.foi, timeStep);
        waveStimulus = DoFreqAnalysis(dataStimulus, data.bands{b}.foi, timeStep);
        [ normalizedMeanDiff restVar stimVar ] = collapsePowerspectrum( waveRest.powspctrm, waveStimulus.powspctrm );
        data.bands{b}.normalizedMeanDiff = normalizedMeanDiff;
        data.bands{b}.restVar = restVar;
        data.bands{b}.stimulusVar = stimVar;
    end
    
    data.electrodeLabels = electrodeLabels;
    data.setDescription = ['VisResp Banded Set ' num2str(x) ' of ' subjectinfo.Name];
    save(wavedataFilename, 'data', '-V7.3');
    
    subjectinfo.Analyses.VisualResponsivenessBanded.Sets(x).WavedataFilename = wavedataFilename;
end

end



%% Various helper methods
function [data, electrodeLabels] = LoadSet(filename)
% Load the data and convert to fieldtrip structure
EEG = pop_loadset( 'filename', filename);
data = eeglab2fieldtrip(EEG, 'preprocessing');
electrodeLabels = data.label;
end

function [wave] = DoFreqAnalysis(data, foi, timeStep)
cfg = [];
cfg.channel    = 'all';
cfg.method     = 'wltconvol';
cfg.output     = 'pow';
cfg.keeptrials = 'no'; % will average all epochs..
cfg.foi        = foi;
cfg.toi        = 1:timeStep:max(data.time{1});
wave = freqanalysis(cfg, data);
end